#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <float.h>
#include <cuda.h>
#include "cublas.h"
#include <time.h>

#define ceil(x)	( x - (x % 1) + 1 )
#define INF 65535

typedef struct graphstruct { // A graph in compressed-adjacency-list (CSR) form
    int nv;            // number of vertices
    int ne;            // number of edges
    int *nbr;          // array of neighbors of all vertices
    int *firstnbr;     // index in nbr[] of first neighbor of each vtx
} graph;


__global__ void baseline_bfs_kernel(const int *,const int *, int *, int , bool *, int );

int read_edge_list (int **tailp, int **headp) {
    int max_edges = 100000000;
    int nedges, nr, t, h;
    *tailp = (int *) calloc(max_edges, sizeof(int));
    *headp = (int *) calloc(max_edges, sizeof(int));
    nedges = 0;
    nr = scanf("%i %i",&t,&h);
    while (nr == 2) {
        if (nedges >= max_edges) {
            printf("Limit of %d edges exceeded.\n",max_edges);
            exit(1);
        }
        (*tailp)[nedges] = t;
        (*headp)[nedges++] = h;
        nr = scanf("%i %i",&t,&h);
    }
    return nedges;
}

graph * graph_from_edge_list (int *tail, int* head, int nedges) {
    graph *G;
    int i, e, v, maxv;
    G = (graph *) calloc(1, sizeof(graph));
    G->ne = nedges;
    maxv = 0;

    // count vertices
    for (e = 0; e < G->ne; e++) {
        if (tail[e] > maxv) maxv = tail[e];
        if (head[e] > maxv) maxv = head[e];
    }
    G->nv = maxv+1;
    G->nbr = (int *) calloc(G->ne, sizeof(int));
    G->firstnbr = (int *) calloc((G->nv)+1, sizeof(int));
    
    // count neighbors of vertex v in firstnbr[v+1],
    for (e = 0; e < G->ne; e++) G->firstnbr[tail[e]+1]++;
    
    // cumulative sum of neighbors gives firstnbr[] values
    for (v = 0; v < G->nv; v++) G->firstnbr[v+1] += G->firstnbr[v];
    
    // pass through edges, slotting each one into the CSR structure
    for (e = 0; e < G->ne; e++) {
        i = G->firstnbr[tail[e]]++;
        G->nbr[i] = head[e];
    }
    // the loop above shifted firstnbr[] left; shift it back right
    for (v = G->nv; v > 0; v--) G->firstnbr[v] = G->firstnbr[v-1];
    G->firstnbr[0] = 0;
    return G;
}


void print_CSR_graph (graph *G) {
    int vlimit = 20;
    int elimit = 50;
    int e,v;
    printf("\nGraph has %d vertices and %d edges.\n",G->nv,G->ne);
    printf("firstnbr =");
    if (G->nv < vlimit) vlimit = G->nv;
    for (v = 0; v <= vlimit; v++) printf(" %d",G->firstnbr[v]);
    if (G->nv > vlimit) printf(" ...");
    printf("\n");
    printf("nbr =");
    if (G->ne < elimit) elimit = G->ne;
    for (e = 0; e < elimit; e++) printf(" %d",G->nbr[e]);
    if (G->ne > elimit) printf(" ...");
    printf("\n\n");
}


int main(int argc, char **argv)
{
	graph *G;
	int *tail, *head ;
	int nedges , nlevels, startvtx ;
	cudaError_t err ;

	//My variables
	int curr ,i, *ctr; 
	bool *finished , *devfinished;
	int *nodes, *edges, *gpulevel , *cpulevel, *testvar; 
	int *testvar_cpu; 

	finished = (bool *) malloc ( 1 * sizeof(bool));
	testvar = (int *) malloc( 1 * sizeof(int));

	nedges = read_edge_list(&tail, &head);
	G = graph_from_edge_list(tail, head, nedges);
	
	cpulevel = (int *)malloc( (G->nv) * sizeof(int));
	cpulevel[0] = 0 ;

	for(i = 1 ; i < G->nv ; i++)
		cpulevel[i] = INF ;

	free(tail);
	free(head);
	print_CSR_graph(G);

	cudaMalloc( (void **)&devfinished, 1 * sizeof(int)) ; 

	curr = 0 ;

	err = cudaMalloc( (void **)&nodes,(G->nv+1) * sizeof(int));
	
	if( err == cudaErrorMemoryAllocation)
	{
		printf("Mem alloc GPU error\n");	
		exit(0);
	}

	err = cudaMalloc( (void **)&edges, (G->ne) * sizeof(int));
	if( err == cudaErrorMemoryAllocation)
	{
		printf("Mem alloc GPU error\n");	
		exit(0);
	}

	err = cudaMalloc( (void **)&gpulevel, (G->nv) * sizeof(int));
	if( err == cudaErrorMemoryAllocation)
	{
		printf("Mem alloc GPU error\n");	
		exit(0);
	}


	dim3 numthreadsPerBlock(16,16) ;
	dim3 numBlocks(ceil((int)sqrt(G->nv)/16),ceil((int)sqrt(G->nv)/16));

	printf("Copying Graph from Host to Device :\n");

	err = cudaMemcpy(nodes,G->firstnbr, (G->nv+1) * sizeof(int) , cudaMemcpyHostToDevice);
	if( err != cudaSuccess)
	{
		printf("Host to Device memory transfer error\n");
		exit(0);
	}

	err = cudaMemcpy(edges,G->nbr, (G->ne) * sizeof(int) ,  cudaMemcpyHostToDevice);
	if( err != cudaSuccess)
	{
		printf("Host to Device memory transfer error\n");
		exit(0);
	}

	err = cudaMemcpy(gpulevel,cpulevel, (G->nv) * sizeof(int) , cudaMemcpyHostToDevice);
	if( err != cudaSuccess)
	{
		printf("Host to Device memory transfer error\n");
		exit(0);
	}

	printf("GPU performing BFS ...\n");

	clock_t start=clock();

	do
	{
	    *finished = true ;

	    cudaMemcpy(devfinished,finished,1 * sizeof(bool), cudaMemcpyHostToDevice);
		if(err != cudaSuccess)
		{
			printf("Host to Device memory transfer error\n");
		}
		else

	    baseline_bfs_kernel<<<numBlocks,numthreadsPerBlock>>>(nodes,edges,gpulevel,curr++,devfinished,G->nv);

	    cudaMemcpy(finished,devfinished,1 * sizeof(bool), cudaMemcpyDeviceToHost);

	}
	while(!(*finished));
 
	printf("Copying result to host ...\n");
	cudaMemcpy(cpulevel , gpulevel, (G->nv) * sizeof(int) , cudaMemcpyDeviceToHost);

	printf("Time elapsed: %f\n", ((double)clock() - start) / CLOCKS_PER_SEC);

    /*if (G->nv < 20) {
	printf("\n  vertex parent  level\n");
	for (v = 0; v < G->nv; v++) printf("%6d%7d%7d\n", v, parent[v], level[v]);
    }*/

	printf("Number of levels = %d\n",curr);

	ctr = (int *)malloc((curr+1) * sizeof(int));

	for(i = 0 ; i < curr ; i++)
		ctr[i] = 0 ; 

	for(i = 0 ; i < (G->nv) ; i++)
	    if(cpulevel[i] < INF)
		ctr[cpulevel[i]]++;

	for(i = 0 ; i < curr ; i++)
		printf("level %d vertices: %d\n",i, ctr[i]);

	return 0 ; 
}

__global__ void baseline_bfs_kernel(const int *nodes,const int *edges, int *levels, int curr, bool *finished, int totalnodes)
{
	int nodeid ,w, i, num_nbrs, row, col;

	col = blockIdx.x * blockDim.x + threadIdx.x ;
        row = blockIdx.y * blockDim.y + threadIdx.y ;	

	nodeid = row * gridDim.x * blockDim.x   + col ;

	if( nodeid < totalnodes)
	{
		if(levels[nodeid] == curr)
		{
			num_nbrs = nodes[nodeid+1] - nodes[nodeid] ;
			
			for( i = 0 ; i < num_nbrs ; i++)
			{
				w = edges[ nodes[nodeid] + i ] ;
				if ( levels[w] == INF)
				{
					*finished = false ;
					levels[w] = curr + 1 ;
				}
			}

		} 
	}
}

